---
title: "Exploring land use prediction errors from area frame survey data"
author: "Supplemental material" 
date: "Raja  Chakir, Thibault  Laurent, Anne  Ruiz-Gazen, Christine  Thomas-Agnan, Céline  Vignes" 
output:
  html_document: 
    toc: true # table of content true
    toc_depth: 2  # upto three depths of headings (specified by #, ## and ###)
    number_sections: true  ## if you want number sections at each table header
 #   theme: united  # many options for theme, this one is my favorite.
    highlight: tango  # specifies the syntax highlighting style
---
<link href="markdown7.css" rel="stylesheet">

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

We provide the data and the R code used in the article "Exploring land use prediction errors from area frame survey data" so that readers may be able to reproduce all the figures, tables and statistics presented in the article with the R software. 

Please cite the use of our resources as: 

Chakir R., Laurent T., Ruiz-Gazen A., Thomas-Agnan C. and Vignes C. (2018). Exploring land use prediction errors from area frame survey data. *CSBIGS*.

# Prerequisites #

Packages needed:
```{r, eval = F}
install.packages(c("rgdal", "maptools", "RColorBrewer", "classInt", "spdep"))
```

Loading packages:
```{r, }
require("rgdal")  # import spatial data
require("maptools") # leglabs function
library("RColorBrewer") # brewer.pal function
require("classInt") # function classIntervals
require("spdep") # spatial econometric analysis
```

Information about the current R session :
```{r}
sessionInfo()
```

# Importing the data #

## Point levels data ## 

The main data (point levels data) is given in the csv file "./data/point_level.csv". It contains 502205 observations and 21 variables:

* **Ui**: value of the categorical variable "simulated land use" (1 = urban, 2 = farming, 3 = forests, 4 = pastures, 5 = natural land),
* **pi1**, **pi2**, **pi3**, **pi4**, **pi5**: the DGP probabilities obtained using the regression tree (see Figure 2),
* **clc2**: corine land cover (in 13 classes, see http://www.eea.europa.eu/publications/COR0-landcover),
* **altitude**: altitude in meters (source: BDAlti, IGN),
* **land.price**: land price in euros/km2 (source: Agreste),
* **pop**: population density in inhabitant/km2 (source: INSEE),
* **INSEE_COM**: the id of the commune,
* **teruti**: a boolean with **TRUE** if the point is used in the French Teruti-Lucas database,
* **numseg**: the id of the segment,
* **A6**, **A5**, **A4**, **A3**, **A2**, **A03**, **A02**, **A01** : the id of the regular square grids at several scales (see Chakir et al (2016), Spatial scale in land use models: Application to the Teruti-Lucas survey, *Spatial Statistics*, 18(A), 246--262, doi:10.1016/j.spasta.2016.06.009).  


For importing the point level data:
```{r}
point.level <- read.csv2("./data/point_level.csv")
```

For transforming the data into Spatial object:
```{r}
coordinates(point.level) <- ~ x + y
proj4string(point.level) <- CRS("+proj=utm +zone=31 +ellps=GRS80 
                                +units=m +no_defs")
```

Please, note that he CRS used for all the geographical coordinates is:
```{r, eval =F}
CRS("+proj=utm +zone=31 +ellps=GRS80 +units=m +no_defs")
```

## The boundaries of different administrative areas ## 

The boundaries of different administrative areas are included in the following folders: 

* **./data/cities**: a shapefile containing the boundaries of the communes of France in Midi-Pyrénées region (source: IGN),
* **./data/dept**: a shapefile containing the boundaries of the departements in Midi-Pyrénées region (source: IGN),
* **./data/segments**: a shapefile containing the boundaries of the segment levels data.

For importing the shapes of the departements:
```{r}
dept <- readOGR(dsn = "./data/dept", layer = "dept")
proj4string(dept)
```

For importing the shapes of the cities:
```{r}
cities <- readOGR(dsn = "./data/cities", layer = "cities")
proj4string(cities)
```

For importing the shapes of the segments
```{r}
segment.level <- readOGR(dsn = "./data/segments", layer = "segments")
proj4string(segment.level)
```

# Preparing the data #

## Preparation of point level data ## 

**Step 1)** dichotomization of **Ui** by creating variables **di1**, **di2**, **di3**, **di4**, **di5**:
```{r}
di <- model.matrix(~ as.factor(point.level@data$Ui) - 1)
point.level@data[,c("di1", "di2", "di3", "di4", "di5")] <- di
```

**Step 2)** we create the absolute response error $|p_{ik} - d_{ik}|$:
```{r}
point.level@data$erreur1.abs <- abs(point.level@data$pi1 - point.level@data$di1)
point.level@data$erreur2.abs <- abs(point.level@data$pi2 - point.level@data$di2)
point.level@data$erreur3.abs <- abs(point.level@data$pi3 - point.level@data$di3)
point.level@data$erreur4.abs <- abs(point.level@data$pi4 - point.level@data$di4)
point.level@data$erreur5.abs <- abs(point.level@data$pi5 - point.level@data$di5)
```

**Step 3)** we create the relative response error $\frac{|p_{ik} - d_{ik}|}{p_{ik}}$:
```{r}
point.level@data$erreur1.rel <- point.level@data$erreur1.abs/point.level@data$pi1
point.level@data$erreur2.rel <- point.level@data$erreur2.abs/point.level@data$pi2
point.level@data$erreur3.rel <- point.level@data$erreur3.abs/point.level@data$pi3
point.level@data$erreur4.rel <- point.level@data$erreur4.abs/point.level@data$pi4
point.level@data$erreur5.rel <- point.level@data$erreur5.abs/point.level@data$pi5
```


**Step 4)** we create the Gini-Simpson impurity index **gi** (noted $gs_i$ in the article):
```{r}
point.level@data$gi <- 1 - rowSums(point.level@data[, c("pi1","pi2","pi3",
                                                        "pi4","pi5")]^2)
```


## Preparation of areal level data ## 

**Step 1)** we create the average of point level Gini-Simpson index **gi** (noted $\overline{gs}_g$ in the article):
```{r}
res <- aggregate(point.level@data$gi, 
                 by = list(numseg = point.level@data$numseg), FUN = mean)
colnames(res)[2] <- "gi"
segment.level <- merge(segment.level, res, by = "numseg")
```

**Step 2)** we create the average at areal level of **pik**: **pk.pA** (noted $\bar{g}_{gk}$ in the article):
```{r}
res <- aggregate(point.level@data[,c("pi1", "pi2", "pi3", "pi4", "pi5")],
                 by = list(numseg = point.level@data$numseg), FUN = mean)
colnames(res)[2:6] <- c("p1.pA", "p2.pA", "p3.pA", "p4.pA", "p5.pA")
segment.level <- merge(segment.level, res, by = "numseg")
```

**Step 3)** we create the average at areal level of **dik**: **pk.dA** (noted $\bar{d}_{gk}$ in the article):
```{r}
res <- aggregate(point.level@data[,c("di1", "di2", "di3", "di4", "di5")],
                 by = list(numseg = point.level@data$numseg), FUN = mean)
colnames(res)[2:6] <- c("p1.dA", "p2.dA", "p3.dA", "p4.dA", "p5.dA")
segment.level <- merge(segment.level, res, by = "numseg")
```

**Step 4)** we create the average point level absolute response error:
```{r}
res <- aggregate(point.level@data[, c("erreur1.abs", "erreur2.abs", 
                                      "erreur3.abs", "erreur4.abs", 
                                      "erreur5.abs")],
                 by = list(numseg = point.level@data$numseg), FUN = mean)
segment.level <- merge(segment.level, res, by = "numseg")
```


**Step 5)** we create the areal level absolute response error:
```{r}
segment.level@data$erreur1.abs.areal <- abs(segment.level@data$p1.dA - 
                                              segment.level@data$p1.pA)
segment.level@data$erreur2.abs.areal <- abs(segment.level@data$p2.dA - 
                                              segment.level@data$p2.pA)
segment.level@data$erreur3.abs.areal <- abs(segment.level@data$p3.dA - 
                                              segment.level@data$p3.pA)
segment.level@data$erreur4.abs.areal <- abs(segment.level@data$p4.dA - 
                                              segment.level@data$p4.pA)
segment.level@data$erreur5.abs.areal <- abs(segment.level@data$p5.dA - 
                                              segment.level@data$p5.pA)
rm(res)
```

# Reproducing the figures, tables and comments #

### Figure 1 ###

```{r, eval = F}
op <- par(mar = c(0.5, 0.5, 0.5, 0.5))
plot(cities, axes = T, xlim = c(370000, 380000), ylim = c(4825000, 4835000),
     border = "grey", xaxt = "n", yaxt = "n")
plot(segment.level, lwd = 2, add = T)
pal <- c("red", "gray")
plot(point.level, col = pal[ifelse(point.level@data$teruti, 1, 2)], cex = 0.3,
     add = T)
par(op)
rm(pal)
```

![](figures/figure1.png){width=50%}

Figure 1: Areal level grid (in black) and points (Teruti-Lucas locations in red) in Toulouse area. 

### Figure 2 ###

The figure 2 has been obtained on confidential data. Hence, codes are not given in here.

![](figures/figure2.png){width=80%}

Figure 2: Classification tree chosen for the DGP.

### Figure 3 ### 

```{r, eval = F}
op <- par(mfrow = c(2, 3), mar = c(2, 4, 1.5, 0.5))
plot(table(round(point.level@data$pi1, 3)), las = 1, ylab = "Frequency", 
     lwd = 5, lend = "square", xlim = c(0, 0.85), cex.axis = 0.8)
title(expression(paste("Urban: ",p[i1])), cex.main = 1.5)
plot(table(round(point.level@data$pi2, 3)), las = 1, ylab = "Frequency", 
     lwd = 5, lend = "square", xlim = c(0, 0.85), cex.axis = 0.8)
title(expression(paste("Farming: ", p[i2])), cex.main = 1.5)
plot(table(round(point.level@data$pi3, 3)), las = 1, ylab = "Frequency", 
     lwd = 5, lend = "square", xlim = c(0, 0.85), cex.axis = 0.8)
title(expression(paste("Forests: ", p[i3])), cex.main = 1.5)
plot(table(round(point.level@data$pi4, 3)), las = 1, ylab = "Frequency", 
     lwd = 5, lend = "square", xlim = c(0, 0.85), cex.axis = 0.8)
title(expression(paste("Pastures: ", p[i4])), cex.main = 1.5)
plot(table(round(point.level@data$pi5, 3)), las = 1, ylab = "Frequency", 
     lwd = 5, lend = "square", xlim = c(0, 0.85), cex.axis = 0.8)
title(expression(paste("Natural land: ", p[i5])), cex.main = 1.5)
par(op)
```

![](figures/figure3.png){width=80%}

Figure 3: Bar charts of $p_{ik}$.


### Comments in the text (in section 3.2. The DGP probabilities) ###

```{r, eval =F}
knots(ecdf(point.level@data$pi1))
ecdf(point.level@data$pi1)(0.072)
```
Comment: ``Most of them are below 0.10 with $95.4\%$ of values less than $0.072$ for $p_{i1}$''.

```{r, eval =F}
knots(ecdf(point.level@data$pi5))
ecdf(point.level@data$pi5)(0.084)
```
Comment: ``$90.0\%$ of the values less than 0.084 for $p_{i5}$''.

```{r, eval =F}
round(table(round(point.level@data$pi2, 3) == 0.726)/
        length(point.level@data$pi2)*100, 1)
```
Comment: ``$23.7\%$ of the values equal to 0.726 for farming''.

```{r, eval =F}
round(table(round(point.level@data$pi3, 3) == 0.822)/
        length(point.level@data$pi3)*100, 1)
```
Comment: ``$26.0\%$ of the values equal to 0.822 for forests''.

```{r, eval =F}
round(table(round(point.level@data$pi4, 3) == 0.570)/
        length(point.level@data$pi4)*100, 1)
```
Comment: $21.9\%$ of the values equal to 0.570 for pastures''.

### Figure 4 ###

```{r, eval =F}
op <- par(mfrow = c(2, 3), mar = c(4, 4, 1.8, 0.5))
plot(ecdf(point.level@data$pi1), main = "", xlab = expression(p[i1]),
     ylab = "ECDF", cex = 0.8, xlim = c(0, 0.85))
title("Urban", cex.main = 1.5)
plot(ecdf(point.level@data$pi2), main = "", xlab = expression(p[i2]),
     ylab = "ECDF", cex = 0.8, xlim = c(0, 0.85))
title("Farming", cex.main = 1.5)
plot(ecdf(point.level@data$pi3), main = "", xlab = expression(p[i3]),
     ylab = "ECDF", cex = 0.8, xlim = c(0, 0.85))
title("Forests", cex.main = 1.5)
plot(ecdf(point.level@data$pi4), main = "", xlab = expression(p[i4]),
     ylab = "ECDF", cex = 0.8, xlim = c(0, 0.85))
title("Pastures", cex.main = 1.5)
plot(ecdf(point.level@data$pi5), main = "", xlab = expression(p[i5]),
     ylab = "ECDF", cex = 0.8, xlim = c(0, 0.85))
title("Natural land", cex.main = 1.5)
par(op)

```

![](figures/figure4.png){width=80%}

Figure 4: Empirical cumulative functions of $p_{ik}$.

### Figure 5 ###


```{r, eval = F}
gris <- brewer.pal(5, "Greys")
bk <- c(0, 0.20, 0.40, 0.60, 0.80, 1)
op <- par(mfrow = c(2, 3), mar = c(0, 0, 2, 0), oma = c(0, 0, 0, 0))
plot(segment.level, border = "gray", 
  col = gris[findInterval(segment.level@data$p1.pA, bk, all.inside = T)])
plot(dept, add = T)
title("Urban", cex.main = 2)
plot(segment.level, border = "gray", 
 col = gris[findInterval(segment.level@data$p2.pA, bk, all.inside = T)])
plot(dept, add = T)
title("Farming", cex.main = 2)
plot(segment.level, border = "gray", 
 col = gris[findInterval(segment.level@data$p3.pA, bk, all.inside = T)])
plot(dept, add = T)
title("Forests", cex.main = 2)
plot(segment.level, border = "gray", 
 col = gris[findInterval(segment.level@data$p4.pA, bk, all.inside = T)])
plot(dept, add = T)
title("Pastures", cex.main = 2)
plot(segment.level, border = "gray", 
 col = gris[findInterval(segment.level@data$p5.pA, bk, all.inside = T)])
plot(dept, add = T)
title("Natural land", cex.main = 2)
plot(rep(1, 5), rev(seq(2, 4, 0.5)), axes = FALSE, xlim = c(0, 3),
 ylim = c(1, 5), pch = 15, cex = 3, col = gris, xlab = "", ylab = "")
text(1.8, 4.5, expression(bar(p)[gk]), cex = 2)
text(rep(1.2, 5), rev(seq(2, 4, 0.5)), 
     c("[0-0.2[", "[0.2-0.4[", "[0.4-0.6[", "[0.6-0.8[", "[0.8-1]"),
 pos = 4, offset = 0.2, cex = 1.5)
par(op)
rm(gris, bk)
```

![](figures/figure5.png){width=80%}

Figure 5: DGP probabilities $\bar{p}_{gk}$.

### Figure 6 ### 

```{r, eval = F}
op <- par(mar = c(4, 4, 0.5, 0.5))
plot(table(round(point.level@data$gi,3)), las = 1, lwd = 5, lend = "square",
 cex.axis = 0.7, ylab = "Frequency",
 xlab = expression(paste("Point level Gini-Simpson impurity index (", gs[i],")")))
par(op)
```

![](figures/figure6.png){width=60%}

Figure 6: Bar chart of point level GS impurity index ($gs_i$).

### Table 2 ### 

**Step 1)** Recoding of **gi** in 8 groups according to the bar chart and the contingency table: **gi.cl** variable
```{r}
point.level@data$gi.cl <- cut(point.level@data$gi,
                              breaks = c(0.3, 0.32, 0.38, 0.46, 0.463, 
                                         0.51, 0.61, 0.62, 0.722))
levels(point.level@data$gi.cl) <- c("0.314", "0.370", "0.450", "0.462", 
                                    "0.509", "0.604", "0.619", "0.682-0.721")
```

```{r, eval = F}
table(point.level@data$gi.cl)
```

**Step 2)** Contingency table landuse (**Ui** vs **Gini**) in 9 groups (**gi.cl**)
```{r, eval = F}
addmargins(table(point.level@data$Ui, point.level@data$gi.cl))
```

Obtaining row percents:
```{r, eval = F}
round(addmargins(prop.table(table(point.level@data$Ui, point.level@data$gi.cl), 
                            margin = 1))*100, 1)
```

Obtaining column percents:
```{r, eval = F}
round(addmargins(prop.table(table(point.level@data$Ui, point.level@data$gi.cl),
                            margin = 2))*100, 1)
```


The column percent table is used to complete Table 2: Characteristics of groups according to the $gs_i$ value (2 first columns).


### Figure 7 ### 

```{r, eval = F}
op <- par(mar = c(4.5, 4.5, 0.5, 0.5))
plot(factor(point.level@data$Ui), point.level@data$gi, xlab = "land use",
 ylim = c(0.3, 0.75), axes = F,
 ylab = expression(paste("Point level Gini-Simpson impurity index   ", 
                         gs[i])))
axis(1, at = seq(1, 5), labels = c("urban", "farming", "forests", 
                                   "pastures", "natural land"))
axis(2)
box()
par(op)
```

![](figures/figure7.png){width=60%}

Figure 7 : Boxplots of $gs_i$ by land use.

### Figure 8 ### 

```{r, eval = F}
blue.pal <- brewer.pal(5, "Blues")
bk <- round(classIntervals(segment.level@data$gi, n = 5, style = "kmeans")$brks,
            digits = 2)
op <- par(mar = c(0, 0, 0, 0))
plot(segment.level, border = "gray",
     col = blue.pal[findInterval(segment.level@data$gi, 
                                 bk, all.inside = T)])
legend("topleft", legend = leglabs(bk), fill = blue.pal, inset = 0.05, 
       title = expression(bar(gs)[g]))
par(op)
```


![](figures/figure8.png){width=60%}

Figure 8: Map of $\bar{gs}_g$, the mean of $gs_i$ at areal level.

### Comments given in Section 4.2. ###

Spatial analysis of the areal level Gini-Simpson index **gi** ($\bar{gs}_g$ in the article).

**Step 1)** Standardization of gi
```{r}
x <- as.vector(scale(segment.level$gi))
```

**Step 2)** Creation of the neighborhood matrix (8 nearest neighbors), normalized
```{r}
colw <- nb2listw(knn2nb(knearneigh(coordinates(segment.level), k = 8)), 
                 style = "W")
```

**Step 3)** Compute the Moran test
```{r}
# 1. "free sampling model" version
# highly significant (p-value < 2.2e-16): presence of spatial autocorrelation
# normalized Moran I = 0.59
moran.test(x, colw, randomisation = FALSE)
# 2. permutation test: "randomization model" version
# significant at 1%: presence of spatial autocorrelation, Moran I = 0.59
set.seed(1234)
moran.mc(x, colw, nsim = 99)

```


### Figure 9 ### 

```{r, eval = F}
op <- par(mar = c(4, 4.5, 0.5, 0.5))
res.moran <- moran.plot(x, colw, xlab = expression(bar(gs)[g]),
                        ylab = expression(paste("W x ", bar(gs)[g])),
                        labels = F, pch = 20, cex = 0.6, cex.axis = 0.8, 
                        cex.lab = 0.9)
par(op)
```

![](figures/figure9.png){width=60%}

Figure 9: Moran scatterplot of the Gini-Simpson index $\bar{gs}_g$ (after standardization).


### Figure 10 ### 

```{r, eval = F}
col.quad <- c("brown3", "bisque", "darkcyan", "darkseagreen1")
wx <- lag.listw(colw, x)

x1 <- which((x > mean(x)) & (wx > mean(wx)))	# HH
x2 <- which((x > mean(x)) & (wx < mean(wx)))	# HL
x3 <- which((x < mean(x)) & (wx < mean(wx)))	# LL
x4 <- which((x < mean(x)) & (wx > mean(wx)))	# LH

i.col <- rep("white", length(x))
i.col[x1] <- col.quad[1]
i.col[x2] <- col.quad[2]
i.col[x3] <- col.quad[3]
i.col[x4] <- col.quad[4]

influ <- rep("", length(x))
influ[apply(res.moran$is.inf, 1, any)] <- 
  row.names(segment.level)[apply(res.moran$is.inf, 1, any)]

op <- par(mar = c(0, 0, 0, 0))
plot(segment.level, col = i.col, cex = c(0.1, 1)[as.integer(influ == "") + 1],
     border = "gray")
plot(dept, lwd = 2, add = T)
legend("topleft", c("H-H", "H-L", "L-L", "L-H"), fill = col.quad, 
       title = "Cluster", inset = 0.05)
par(op)
```

![](figures/figure10.png){width=60%}

Figure 10: Clusters of the Gini-Simpson index $\bar{gs}_g$.

### Figure 11 ### 

```{r, eval = F}
resI <- localmoran(x, colw)
resI1 <- which((resI[, 5] < 0.05) & (x > mean(x)) & (wx > mean(wx)))	# HH
resI2 <- which((resI[, 5] < 0.05) & (x > mean(x)) & (wx < mean(wx)))	# HL
resI3 <- which((resI[, 5] < 0.05) & (x < mean(x)) & (wx < mean(wx)))	# LL
resI4 <- which((resI[, 5] < 0.05) & (x < mean(x)) & (wx > mean(wx)))	# LH

p.col <- rep("white", length(x))
p.col[resI1] <- col.quad[1]
p.col[resI2] <- col.quad[2]
p.col[resI3] <- col.quad[3]
p.col[resI4] <- col.quad[4]

op <- par(mar = c(0, 0, 0, 0))
plot(segment.level, col = p.col, border = "gray")
plot(dept, lwd = 2, add = T)
legend("topleft", c("H-H", "L-L"), fill = c("brown3", "darkcyan"), 
       title = "Cluster", inset = 0.05)
par(op)

rm(x, colw, wx, res.moran, col.quad, x1, x2, x3, x4, i.col, influ, resI, resI1,
   resI2, resI3, resI4, p.col)
```

![](figures/figure11.png){width=60%}

Figure 11: Segments with significant LISA for the $\bar{gs}_g$.

### Table 3 ### 

**Step 1)** Absolute error when $d_{ik}=0$ (1st part of the table)
```{r, eval = F}
round(summary(point.level@data$erreur1.abs[point.level@data$di1 == 0]), 2)
round(summary(point.level@data$erreur2.abs[point.level@data$di2 == 0]), 2)
round(summary(point.level@data$erreur3.abs[point.level@data$di3 == 0]), 2)
round(summary(point.level@data$erreur4.abs[point.level@data$di4 == 0]), 2)
round(summary(point.level@data$erreur5.abs[point.level@data$di5 == 0]), 2)

length(point.level@data$erreur1.abs[point.level@data$di1 == 0])
length(point.level@data$erreur2.abs[point.level@data$di2 == 0])
length(point.level@data$erreur3.abs[point.level@data$di3 == 0])
length(point.level@data$erreur4.abs[point.level@data$di4 == 0])
length(point.level@data$erreur5.abs[point.level@data$di5 == 0])
```


**Step 2)** Absolute error when $d_{ik}=1$ (2nd part of the table)
```{r, eval = F}
round(summary(point.level@data$erreur1.abs[point.level@data$di1 == 1]), 2)
round(summary(point.level@data$erreur2.abs[point.level@data$di2 == 1]), 2)
round(summary(point.level@data$erreur3.abs[point.level@data$di3 == 1]), 2)
round(summary(point.level@data$erreur4.abs[point.level@data$di4 == 1]), 2)
round(summary(point.level@data$erreur5.abs[point.level@data$di5 == 1]), 2)

length(point.level@data$erreur1.abs[point.level@data$di1==1])
length(point.level@data$erreur2.abs[point.level@data$di2==1])
length(point.level@data$erreur3.abs[point.level@data$di3==1])
length(point.level@data$erreur4.abs[point.level@data$di4==1])
length(point.level@data$erreur5.abs[point.level@data$di5==1])
```

**Step 3)** Relative error when $d_{ik}=1$ (3rd part of the table)

```{r, eval = F}
round(summary(point.level@data$erreur1.rel[point.level@data$di1 == 1]), 2)
round(summary(point.level@data$erreur2.rel[point.level@data$di2 == 1]), 2)
round(summary(point.level@data$erreur3.rel[point.level@data$di3 == 1]), 2)
round(summary(point.level@data$erreur4.rel[point.level@data$di4 == 1]), 2)
round(summary(point.level@data$erreur5.rel[point.level@data$di5 == 1]), 2)

length(point.level@data$erreur1.rel[point.level@data$di1 == 1])
length(point.level@data$erreur2.rel[point.level@data$di2 == 1])
length(point.level@data$erreur3.rel[point.level@data$di3 == 1])
length(point.level@data$erreur4.rel[point.level@data$di4 == 1])
length(point.level@data$erreur5.rel[point.level@data$di5 == 1])
```

Table 3: Descriptives statistics of absolute and relative error at point level.

### Figure NOT included in this article but in Chakir et al (2016c) ### 

``Another plot for exploring the absolute response error is proposed by Haaf et al. (2014) and can be found in Chakir et al.
(2016b) for our data set of interest. It is called the Cumulative Distribution Function of Error Tolerance (CDFET) and consists in the empirical cumulative distribution function of the response error for each land use.''

```{r, eval = F}
CDFET1 <- ecdf(point.level@data$erreur1.abs)
CDFET2 <- ecdf(point.level@data$erreur2.abs)
CDFET3 <- ecdf(point.level@data$erreur3.abs)
CDFET4 <- ecdf(point.level@data$erreur4.abs)
CDFET5 <- ecdf(point.level@data$erreur5.abs)

op <- par(mar = c(4, 4, 0.5, 0.5))
plot(CDFET1, verticals = T, col = "red", main = "", xlim = c(0, 1),
     ylab = "CDFET", cex = 0.8, 
     xlab = expression(paste("Pointwise absolute response error: ",
                             group("|",d[ik] - p[ik], "|"))))
plot(CDFET2, verticals = T, col = "yellow", main = "", cex = 0.8, add = T)
plot(CDFET3, verticals = T, col = "darkgreen", main = "", cex = 0.8, add = T)
plot(CDFET4, verticals = T, col = "palegreen", main = "", cex = 0.8, add = T)
plot(CDFET5, verticals = T, col = "grey", main = "", cex = 0.8, add = T)
legend("right", legend = c("urban", "farming", "forests", "pastures", 
                           "natural land"), lty = 1, 
       col = c("red", "yellow", "darkgreen", "palegreen", "grey"), 
       inset = 0.01)
par(op)

rm(CDFET1, CDFET2, CDFET3, CDFET4, CDFET5)
```

![](figures/figurebis.png){width=90%}

### Figure 12 ### 

Figure 12 : Absolute response error vs Gini-Simpson index $g_i$.

1st row of the figure:

```{r, eval = F}
op <- par(mfrow = c(1, 5), oma = c(0, 0.3, 1.5, 0), mar = c(1.1, 2.3, 3.5, 0),
          mai = c(0.05, 0.13, 0.3, 0.0), omi = c(0, 0.5, 0.05, 0))
bp <- with(subset(point.level@data, di1 == 1), 
           boxplot(erreur1.abs ~ gi.cl, xaxt = "n", las = 1, 
                   ylim = c(0, 1.1), cex.axis = 1.5))
mtext("urban observed", side = 2, line = 3, cex = 1.3, font = 2)
text(2.5, 1.1, expression(paste(group("|", d[i1] - p[i1], "|"), " when ",
                                d[i1], "=1")), cex = 1.5)
axis(1, at = 1:8, labels = FALSE)
title("error on urban", cex.main = 1.9)
bp <- with(subset(point.level@data, di1 == 1), 
           boxplot(erreur2.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i2] - p[i2], "|"), " when ",
                                d[i1],"=1")), cex = 1.5)
title("error on farming", cex.main = 1.9)
bp <- with(subset(point.level@data, di1 == 1), 
           boxplot(erreur3.abs ~ gi.cl, xaxt = "n", yaxt = "n",
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i3] - p[i3], "|"), " when ",
                                d[i1], "=1")), cex = 1.5)
title("error on forests", cex.main = 1.9)
bp <- with(subset(point.level@data, di1 == 1), 
           boxplot(erreur4.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i4] - p[i4], "|"), " when ",
                                d[i1], "=1")), cex = 1.5)
title("error on pastures", cex.main = 1.9)
bp <- with(subset(point.level@data, di1 == 1), boxplot(erreur5.abs ~ gi.cl,
 xaxt = "n", yaxt = "n", las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i5] - p[i5],"|"), " when ",
                                d[i1], "=1")), cex = 1.5)
title("error on natural land", cex.main = 1.9)
par(op)
```

![](figures/figure12a.png){width=90%}

2nd row of the figure:

```{r, eval = F}
op <- par(mfrow = c(1, 5), oma = c(0, 0, 0, 0), mar = c(1.1, 2.0, 1.8, 0),
          mai = c(0.05, 0.1, 0.05, 0.0), omi = c(0, 0.5, 0, 0))
bp <- with(subset(point.level@data, di2 == 1), 
           boxplot(erreur1.abs ~ gi.cl, xaxt = "n", las = 1, 
                   ylim = c(0, 1.1), cex.axis = 1.5))
mtext("farming observed", side = 2, line = 3, cex = 1.3, font = 2)
text(2.5, 1.1, expression(paste(group("|", d[i1] - p[i1], "|"), " when ",
                                d[i2], "=1")), cex = 1.5)
axis(1, at = 1:8, labels = FALSE)
bp <- with(subset(point.level@data, di2 == 1), 
           boxplot(erreur2.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1), cex.axis = 1.5))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i2] - p[i2],"|"), " when ",
                                d[i2], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di2 == 1), 
           boxplot(erreur3.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i3] - p[i3], "|"), " when ",
                                d[i2], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di2 == 1), 
           boxplot(erreur4.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i4] - p[i4], "|"), " when ",
                                d[i2], "=1")), cex = 1.5)
bp <- with(subset(point.level@data,di2==1), 
           boxplot(erreur5.abs ~ gi.cl, xaxt = "n", 
                   yaxt = "n", las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i5] - p[i5], "|"), " when ",
                                d[i2], "=1")), cex = 1.5)
par(op)
```

![](figures/figure12b.png){width=90%}

3rd row of the figure:

```{r, eval = F}
op <- par(mfrow = c(1, 5), oma = c(0, 0, 0, 0), mar = c(1.1, 2.0, 1.8, 0),
          mai = c(0.05, 0.1, 0.05, 0.0), omi = c(0, 0.5, 0, 0))
bp <- with(subset(point.level@data, di3 == 1), 
           boxplot(erreur1.abs ~ gi.cl, xaxt = "n", las = 1,
                   ylim = c(0, 1.1), cex.axis = 1.5))
mtext("forests observed", side = 2, line = 3, cex = 1.3, font = 2)
text(2.5, 1.1, expression(paste(group("|", d[i1] - p[i1], "|"), " when ",
                                d[i3], "=1")), cex = 1.5)
axis(1, at = 1:8, labels = FALSE)
bp <- with(subset(point.level@data, di3 == 1), 
           boxplot(erreur2.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1), cex.axis = 1.5))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i2] - p[i2],"|"), " when ",
                                d[i3], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di3 == 1), boxplot(erreur3.abs~gi.cl,
 xaxt = "n", yaxt = "n", las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i3] - p[i3],"|"), " when ",
                                d[i1], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di3 == 1), 
           boxplot(erreur4.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i4] - p[i4],"|"), " when ",
                                d[i3], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di3 == 1), 
           boxplot(erreur5.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i5] - p[i5], "|"), " when ",
                                d[i3], "=1")), cex = 1.5)
par(op)
```

![](figures/figure12c.png){width=90%}

4th row of the figure:

```{r, eval = F}
op <- par(mfrow = c(1,5), oma = c(0, 0, 0, 0), mar = c(1.1, 2.0, 1.8, 0),
          mai = c(0.05, 0.1, 0.05, 0.0), omi = c(0, 0.5, 0, 0))
bp <- with(subset(point.level@data, di4 == 1), 
           boxplot(erreur1.abs ~ gi.cl, xaxt = "n", las = 1, 
                   ylim = c(0, 1.1), cex.axis = 1.5))
mtext("pastures observed", side = 2, line = 3, cex = 1.3, font = 2)
text(2.5, 1.1, expression(paste(group("|", d[i1] - p[i1], "|"), " when ",
                                d[i4], "=1")), cex = 1.5)
axis(1, at = 1:8, labels = FALSE)
bp <- with(subset(point.level@data, di4 == 1), 
           boxplot(erreur2.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i2] - p[i2],"|"), " when ",
                                d[i4], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di4 == 1), 
           boxplot(erreur3.abs ~ gi.cl, xaxt = "n", 
                   yaxt = "n", las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i3] - p[i3], "|"), " when ",
                                d[i4], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di4 == 1), 
           boxplot(erreur4.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i4] - p[i4],"|"), " when ",
                                d[i4], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di4 == 1), 
           boxplot(erreur5.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:8, labels = FALSE)
text(2.5, 1.1, expression(paste(group("|", d[i5] - p[i5],"|"), " when ",
                                d[i4], "=1")), cex = 1.5)
par(op)
```

![](figures/figure12d.png){width=90%}

5th row of the figure:

```{r, eval = F}
op <- par(mfrow = c(1, 5), oma = c(0, 0, 0, 0), mar = c(1.8, 2.0, 1.8, 0),
          mai = c(0.7, 0.1, 0.05, 0.0), omi = c(0.15, 0.5, 0, 0))
bp <- with(subset(point.level@data, di5 == 1), 
           boxplot(erreur1.abs ~ gi.cl, xaxt = "n", las = 1, 
                   ylim = c(0, 1.1), cex.axis = 1.5))
mtext("natural land observed", side = 2, line = 3, cex = 1.3, font = 2)
text(1:9, -0.09, labels = levels(point.level@data$gi.cl), srt = 30, adj = 1,
     xpd = TRUE, cex = 1.7)
text(2.5, 1.1, expression(paste(group("|", d[i1] - p[i1], "|"), " when ",
                                d[i5], "=1")), cex = 1.5)
axis(1, at = 1:9, labels = FALSE)
bp <- with(subset(point.level@data, di5 == 1), 
           boxplot(erreur2.abs ~ gi.cl, xaxt = "n", 
                   yaxt = "n", las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:9, labels = FALSE)
text(1:9, -0.09, labels = levels(point.level@data$gi.cl), srt = 30, adj= 1,
     xpd = TRUE, cex=1.7)
text(2.5, 1.1, expression(paste(group("|", d[i2] - p[i2],"|"), " when ",
                                d[i5], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di5 == 1), 
           boxplot(erreur3.abs ~ gi.cl, xaxt = "n", yaxt = "n",
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:9, labels = FALSE)
text(1:9, -0.09, labels = levels(point.level@data$gi.cl),
     srt = 30, adj= 1, xpd = TRUE, cex=1.7)
text(2.5, 1.1, expression(paste(group("|", d[i3] - p[i3], "|"), " when ",
                                d[i5], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di5 == 1), 
           boxplot(erreur4.abs ~ gi.cl, xaxt = "n", yaxt = "n",
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:9, labels = FALSE)
text(1:9, -0.09, labels = levels(point.level@data$gi.cl), srt = 30, adj= 1,
     xpd = TRUE, cex=1.7)
text(2.5, 1.1, expression(paste(group("|", d[i4] - p[i4], "|"), " when ",
                                d[i5], "=1")), cex = 1.5)
bp <- with(subset(point.level@data, di5 == 1), 
           boxplot(erreur5.abs ~ gi.cl, xaxt = "n", yaxt = "n", 
                   las = 1, ylim = c(0, 1.1)))
axis(2, labels = FALSE)
axis(1, at = 1:9, labels = FALSE)
text(1:9, -0.09, labels = levels(point.level@data$gi.cl), srt = 30, adj= 1,
     xpd = TRUE, cex=1.7)
text(2.5, 1.1, expression(paste(group("|", d[i5] - p[i5],"|"), " when ",
                                d[i5], "=1")), cex = 1.5)
par(op)
rm(bp)
```

![](figures/figure12e.png){width=90%}

### Table 4 ### 

Table 4: Descriptive statistics of the two types of absolute response error at areal level

**Step 1)** Average point level absolute response error:
```{r, eval = F}
round(summary(segment.level@data$erreur1.abs), 2)
round(summary(segment.level@data$erreur2.abs), 2)
round(summary(segment.level@data$erreur3.abs), 2)
round(summary(segment.level@data$erreur4.abs), 2)
round(summary(segment.level@data$erreur5.abs), 2)
```


**Step 2)** Aggregated absolute response error:
```{r, eval = F}
round(summary(segment.level@data$erreur1.abs.areal), 2)
round(summary(segment.level@data$erreur2.abs.areal), 2)
round(summary(segment.level@data$erreur3.abs.areal), 2)
round(summary(segment.level@data$erreur4.abs.areal), 2)
round(summary(segment.level@data$erreur5.abs.areal), 2)
```


### Figure 13 ###

```{r, eval = F}
bk <- c(0, 0.1, 0.2, 0.3, 0.4, 0.5)
op <- par(mfrow = c(2, 3), mar = c(0, 0, 3, 0))
plot(segment.level, border = "gray",
     col = blue.pal[findInterval(segment.level@data$erreur1.abs, 
                                 bk, all.inside = T)])
title(expression(paste("Urban: ", sum(group("|", d[i1]-p[i1],"|"), i %in% Gg),
                       " /#", Gg)), cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur2.abs,
                                 bk, all.inside = T)])
title(expression(paste("Farming: ", sum(group("|", d[i2]-p[i2],"|"), i %in% Gg),
                       " /#", Gg)), cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur3.abs,
                                 bk, all.inside = T)])
title(expression(paste("Forests: ", sum(group("|",d[i3]-p[i3],"|"), i %in% Gg),
                       " /#", Gg)), cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur4.abs,
                                 bk, all.inside = T)])
title(expression(paste("Pastures: ", sum(group("|",d[i4]-p[i4],"|"), i %in% Gg),
                       " /#", Gg)), cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur5.abs,
                                 bk, all.inside = T)])
title(expression(paste("Natural land: ", sum(group("|", d[i5]-p[i5],"|"),
                                             i %in% Gg), " /#", Gg)), 
      cex.main = 2)
plot(rep(1, 5), rev(seq(2, 4, 0.5)), axes = FALSE, xlim = c(0, 3), ylim = c(1,5),
     pch = 15, cex = 3, col = blue.pal, xlab = "", ylab = "")
text(1.8, 4.5, "Average point level absolute response error", cex = 1.5)
text(rep(1.2, 5), rev(seq(2, 4, 0.5)), c("[0-0.1[", "[0.1-0.2[", "[0.2-0.3[",
                                         "[0.3-0.4[", "[0.4-0.5]"), 
     pos = 4, offset = 0.2, cex = 1.5)
par(op)
rm(bk)
```

![](figures/figure13.png){width=90%}

Figure 13: Average point level absolute response error.

### Figure 14 ### 

```{r, eval = F}
b1 <- round(classIntervals(segment.level@data$erreur1.abs.areal, n = 5,
                           style = "kmeans")$brks, digits = 2)
b2 <- round(classIntervals(segment.level@data$erreur2.abs.areal, n = 5,
                           style = "kmeans")$brks, digits = 2)
b3 <- round(classIntervals(segment.level@data$erreur3.abs.areal, n = 5,
                           style = "kmeans")$brks, digits = 2)
b4 <- round(classIntervals(segment.level@data$erreur4.abs.areal, n = 5,
                           style = "kmeans")$brks, digits = 2)
b5 <- round(classIntervals(segment.level@data$erreur5.abs.areal, n = 5,
                           style = "kmeans")$brks, digits = 2)
bk <- 1/5*(b1 + b2 + b3 + b4 + b5)
bk[6] <- max(b1[6], b2[6], b3[6], b4[6], b5[6])

op <- par(mfrow = c(2, 3), mar = c(0, 0, 3, 0))
plot(segment.level, border = "gray",
     col = blue.pal[findInterval(segment.level@data$erreur1.abs.areal,
                                 bk, all.inside = T)])
title(expression(paste("Urban: ", group("|", bar(d)[g1] - bar(p)[g1], "|"))),
      cex.main = 2)
plot(segment.level, border = "gray",
     col = blue.pal[findInterval(segment.level@data$erreur2.abs.areal, 
                                 bk, all.inside = T)])
title(expression(paste("Farming: ", group("|", bar(d)[g2] - bar(p)[g2], "|"))),
      cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur3.abs.areal,
                                 bk, all.inside = T)])
title(expression(paste("Forests: ", group("|", bar(d)[g3] - bar(p)[g3], "|"))),
      cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur4.abs.areal, 
                                 bk, all.inside = T)])
title(expression(paste("Pastures: ", group("|", bar(d)[g4] - bar(p)[g4], "|"))),
      cex.main = 2)
plot(segment.level, border = "gray", 
     col = blue.pal[findInterval(segment.level@data$erreur5.abs.areal,
                                 bk, all.inside = T)])
title(expression(paste("Natural land: ", group("|", bar(d)[g5]-bar(p)[g5], "|"))),
      cex.main = 2)
plot(rep(1, 5), rev(seq(2, 4, 0.5)), axes = FALSE, xlim = c(0, 3), ylim = c(1, 5),
     pch = 15, cex = 3, col = blue.pal, xlab = "", ylab = "")
text(1.8, 4.5, "Areal level absolute response error", cex = 1.5)
text(rep(1.2, 5), rev(seq(2, 4, 0.5)), c("[0-0.01[", "[0.01-0.18[",
                                         "[0.18-0.032[", "[0.032-0.05[", 
                                         "[0.05-0.15]"), 
     pos = 4, offset = 0.2, cex = 1.5)
par(op)
rm(b1, b2, b3, b4, b5, bk)
```

![](figures/figure14.png){width=90%}

Figure 14: Areal level absolute response error.

### Figure NOT included in this article but in Chakir et al (2016c) ###

``In Chakir  et  al. (2016b),  the  CDFET  for  the  areal  level  absolute  response  error  is  compared  with
the  CDFET  for  the  average  point  level  response  error.  The  comparison  confirms  once  more  that
the  areal  level  response  error  is  very  small  and  for  all  land  uses.''

```{r, eval = F}
CDFET1_agg <- ecdf(segment.level@data$erreur1.abs)
CDFET2_agg <- ecdf(segment.level@data$erreur2.abs)
CDFET3_agg <- ecdf(segment.level@data$erreur3.abs)
CDFET4_agg <- ecdf(segment.level@data$erreur4.abs)
CDFET5_agg <- ecdf(segment.level@data$erreur5.abs)
CDFET1 <- ecdf(segment.level@data$erreur1.abs.areal)
CDFET2 <- ecdf(segment.level@data$erreur2.abs.areal)
CDFET3 <- ecdf(segment.level@data$erreur3.abs.areal)
CDFET4 <- ecdf(segment.level@data$erreur4.abs.areal)
CDFET5 <- ecdf(segment.level@data$erreur5.abs.areal)

op <- par(mfrow = c(1, 2), mar = c(5, 5, 0.5, 0.5))
plot(CDFET1_agg, verticals = T, col = "red", main = "", cex = 0.8,
     xlab = expression(paste("Average point level absolute response error: 1/#", 
                             Gg, sum(group("|",d[ik]-p[ik],"|"), i %in% Gg))), 
     xlim = c(0, 0.5), ylab = "CDFET")
plot(CDFET2_agg, verticals = T, col = "yellow", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET3_agg, verticals = T, col = "darkgreen", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET4_agg, verticals = T, col = "palegreen", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET5_agg, verticals = T, col = "grey", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET1, verticals = T, col = "red", main = "", cex = 0.8,
     xlab = expression(paste("Areal level absolute response error: ", 
                             group("|", paste(bar(d)[gk]-bar(p)[gk]), "|"))), 
     xlim = c(0, 0.5), ylab = "CDFET")
plot(CDFET2, verticals = T, col = "yellow", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET3, verticals = T, col = "darkgreen", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET4, verticals = T, col = "palegreen", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
plot(CDFET5, verticals = T, col = "grey", main = "", xlim = c(0, 0.5),
     cex = 0.8, add = T)
legend("right", legend = c("urban", "farming", "forests", "pastures", 
                           "natural land"), lty = 1, 
       col = c("red", "yellow", "darkgreen", "palegreen", "grey"))
par(op)

rm(CDFET1, CDFET2, CDFET3, CDFET4, CDFET5, CDFET1_agg, CDFET2_agg, CDFET3_agg,
 CDFET4_agg, CDFET5_agg)
```

![](figures/figurebis2.png){width=85%}

### Figure 15 ###

```{r, eval = F}
op <- par(mfrow = c(2, 3), mar = c(5, 5, 3, 0.5))
plot(segment.level@data$gi, segment.level@data$erreur1.abs, xlim = c(0.3, 0.7),
     ylim = c(0, 0.5), xlab = expression(paste("Aggregated impurity index ", 
                                               bar(gs)[g])),
     ylab = expression(paste("1/#", Gg, sum(group("|", d[i1]-p[i1], "|"), 
                                            i %in% Gg))))
panel.smooth(segment.level@data$gi, segment.level@data$erreur1.abs,
             col.smooth = "blue") 
title(expression(paste("Urban: 1/#", Gg, sum(group("|", d[i1]-p[i1],"|"), 
                                             i %in% Gg))),
      cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur2.abs, xlim = c(0.3, 0.7),
     ylim = c(0, 0.5), xlab = expression(paste("Aggregated impurity index ", 
                                               bar(gs)[g])),
     ylab = expression(paste("1/#", Gg, sum(group("|", d[i2]-p[i2], "|"), 
                                            i %in% Gg))))
panel.smooth(segment.level@data$gi, segment.level@data$erreur2.abs,
             col.smooth = "blue")
title(expression(paste("Farming: 1/#", Gg, sum(group("|", d[i2]-p[i2], "|"),
                                               i %in% Gg))), cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur3.abs, xlim = c(0.3, 0.7),
     ylim = c(0, 0.5), xlab = expression(paste("Aggregated impurity index ", 
                                               bar(gs)[g])),
     ylab = expression(paste("1/#", Gg, sum(group("|", d[i3]-p[i3],"|"), 
                                            i %in% Gg))))
panel.smooth(segment.level@data$gi, segment.level@data$erreur3.abs,
             col.smooth = "blue")
title(expression(paste("Forests: 1/#", Gg, sum(group("|", d[i3]-p[i3], "|"),
                                               i %in% Gg))), cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur4.abs, xlim = c(0.3, 0.7),
     ylim = c(0, 0.5), xlab = expression(paste("Aggregated impurity index ", 
                                               bar(gs)[g])),
     ylab = expression(paste("1/#", Gg, sum(group("|", d[i4]-p[i4],"|"), 
                                            i %in% Gg))))
panel.smooth(segment.level@data$gi, segment.level@data$erreur4.abs,
             col.smooth = "blue")
title(expression(paste("Pastures: 1/#", Gg, sum(group("|", d[i4]-p[i4], "|"),
 i %in% Gg))), cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur5.abs, xlim = c(0.3, 0.7),
     ylim = c(0, 0.5), xlab = expression(paste("Aggregated impurity index ", 
                                               bar(gs)[g])),
     ylab = expression(paste("1/#", Gg, sum(group("|", d[i5]-p[i5], "|"), 
                                            i %in% Gg))))
panel.smooth(segment.level@data$gi, segment.level@data$erreur5.abs,
             col.smooth = "blue")
title(expression(paste("Natural land: 1/#", Gg, sum(group("|", d[i5]-p[i5], "|"),
                                                    i %in% Gg))), cex.main = 1.5)
par(op)
```

![](figures/figure15.png){width=90%}

Figure 15: Average point level absolute response error vs aggregated Gini-Simpson impurity index $\bar{gs}_g$.


### Figure 16 ###

```{r, eval = F}
op <- par(mfrow = c(2, 3), mar = c(5, 5, 3, 0.5))
plot(segment.level@data$gi, segment.level@data$erreur1.abs.areal, 
     xlim = c(0.3, 0.7), ylim = c(0, 0.15), 
     xlab = expression(paste("Aggregated impurity index ", bar(gs)[g])),
     ylab = expression(group("|", bar(d)[g1]-bar(p)[g1], "|")))
panel.smooth(segment.level@data$gi, segment.level@data$erreur1.abs.areal,
             col.smooth = "blue") 
title(expression(paste("Urban: ", group("|", bar(d)[g1]-bar(p)[g1], "|"))),
      cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur2.abs.areal, 
     xlim = c(0.3, 0.7), ylim = c(0, 0.15), 
     xlab = expression(paste("Aggregated impurity index ", bar(gs)[g])),
     ylab = expression(group("|", bar(d)[g2]-bar(p)[g2], "|")))
panel.smooth(segment.level@data$gi, segment.level@data$erreur2.abs.areal,
             col.smooth = "blue")
title(expression(paste("Farming: ", group("|", bar(d)[g2]-bar(p)[g2], "|"))),
      cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur3.abs.areal, 
     xlim = c(0.3, 0.7), ylim = c(0, 0.15), 
     xlab = expression(paste("Aggregated impurity index ", bar(gs)[g])),
     ylab = expression(group("|", bar(d)[g3]-bar(p)[g3], "|")))
panel.smooth(segment.level@data$gi, segment.level@data$erreur3.abs.areal,
             col.smooth = "blue")
title(expression(paste("Forests: ", group("|", bar(d)[g3]-bar(p)[g3], "|"))),
      cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur4.abs.areal, 
     xlim = c(0.3, 0.7), ylim = c(0, 0.15), 
     xlab = expression(paste("Aggregated impurity index ", bar(gs)[g])),
     ylab = expression(group("|", bar(d)[g4]-bar(p)[g4], "|")))
panel.smooth(segment.level@data$gi, segment.level@data$erreur4.abs.areal,
             col.smooth = "blue")
title(expression(paste("Pastures: ", group("|", bar(d)[g4]-bar(p)[g4], "|"))),
      cex.main = 1.5)
plot(segment.level@data$gi, segment.level@data$erreur5.abs.areal, 
     xlim = c(0.3, 0.7), ylim = c(0, 0.15), 
     xlab = expression(paste("Aggregated impurity index ", bar(gs)[g])),
     ylab = expression(group("|", bar(d)[g5]-bar(p)[g5], "|")))
panel.smooth(segment.level@data$gi,segment.level@data$erreur5.abs.areal,
             col.smooth = "blue")
title(expression(paste("Natural land: ", 
                       group("|", bar(d)[g5]-bar(p)[g5], "|"))),
      cex.main = 1.5)
par(op)
```

![](figures/figure16.png){width=90%}

Figure 16: Areal absolute response error vs aggregated Gini-Simpson impurity index $\bar{gs}_g$.

### Annex ###

Simulation by multinomial random draw with probas **pi1**, **pi2**, **pi3**, **pi4**, **pi5**
and creation of new variables: **di1**, **di2**, **di3**, **di4**, **di5** (and then **Ui**).

**Step 1)** we identify the leafs of the tree :
```{r}
leaf.tree <- unique(point.level@data[, c("pi1", "pi2", "pi3", "pi4", "pi5")])
```

**Step 2)** for each point, we give a number of leaf between 1 and the number of leaf :
```{r}
res.tree <-  vector("integer", nrow(point.level))
for(k in 1:nrow(leaf.tree))
  res.tree[point.level@data[, "pi1"] == leaf.tree[k, "pi1"]] <- k
```


**Step 3)** simulation :
```{r}
res.sim <- matrix(0, nrow(point.level), 5)
for(k in 1:nrow(leaf.tree))
  res.sim[res.tree == k, ] <- 
  t(rmultinom(length(which(point.level@data[, "pi1"] == leaf.tree$pi1[k])),
              size = 1, prob = leaf.tree[k, c("pi1", "pi2", "pi3", "pi4", "pi5")]))
```

**Step 4)** we obtain **Ui**:
```{r}
Ui <- apply(res.sim, 1, which.max)
```